A computer program listing of the preferred embodiment of the present invention has been separately submitted on a compact disc entitled xe2x80x9cComputer Program Listing Appendix.xe2x80x9d The computer program listing contained therein is incorporate by reference herein in its entirety.
1. Field of the Invention
The invention relates to the exploration of hydrocarbon wells. More particularly, the invention relates to methods for processing waveforms detected by a downhole tool, particularly a sonic logging tool.
2. State of the Art
Sonic logging is a well developed art, and details relating to sonic logging tools and techniques are set forth in xe2x80x9cPorosity Logsxe2x80x9d; Schlumberger Log Interpretation Principles/Applications, Chapter 5, Schlumberger Educational Services, Texas (1987); A. Kurkjian, et al., xe2x80x9cSlowness Estimation from Sonic Logging Waveformsxe2x80x9d, Geoexploration, Vol. 277, pp.215-256 (1991); C. F. Morris et al., xe2x80x9cA New Sonic Array Tool for Full Waveform Logging,xe2x80x9d SPE-13285, Society of Petroleum Engineers (1984); A. R. Harrison et al., xe2x80x9cAcquisition and Analysis of Sonic Waveforms From a Borehole Monopole and Dipole Source . . . xe2x80x9d SPE 20557, pp. 267-282. (Set. 1990); C. V. Kimball and T. L. Marzetta, xe2x80x9cSemblance Processing of Borehole Acoustic Array Dataxe2x80x9d, Geophysics, Vol. 49, pp. 274-281 (March 1984); U.S. Pat. No. 4,131,875 to Ingram; and U.S. Pat. No. 4,594,691 to Kimball et al., all of which are hereby incorporated by reference herein in their entireties. A sonic logging tool typically includes a sonic source (transmitter), and a plurality of receivers which are spaced apart by several inches or feet. In the borehole arts, a sonic signal is transmitted from a sonic source and received at the receivers of the borehole tool which are spaced apart from the sonic source, and measurements are made every few inches as the tool is drawn up the borehole. The sonic signal from the transmitter or source enters the formation adjacent the borehole, and the arrival times and perhaps other characteristics of the receiver responses are recorded. Typically, compressional (P-wave), shear (S-wave) and Stoneley (fluid) arrivals and waves are detected by the receivers and are processed either downhole or uphole. The information which is recorded is typically used to find formation parameters such as formation slowness (the inverse of sonic speed) and semblance, from which pore pressure, porosity, and other determinations can be made. In certain tools such as the DSI (Dipole Sonic Imager) tool (a trademark of Schlumberger), the sonic signals may even be used to image the formation.
Many different techniques for processing the sonic wave signals are known in the art in order to obtain information regarding the borehole and/or formation. Typically, the processing involves digitizing the received signal at a desired sampling rate and then processing the digitized samples according to desired techniques. Examples may be found in U.S. Pat. No. 4,594,691 to Kimball et al. and the references cited therein, as well as in articles such as A. R. Harrison et al., xe2x80x9cAcquisition and Analysis of Sonic Waveforms From a Borehole Monopole and Dipole Source . . . xe2x80x9d SPE 20557, pp. 267-282 (Sept. 1990).
For some time now, compressional slowness has been computed using Slowness-Time Coherence (STC) processing. C. V. Kimball and T. L. Marzetta, xe2x80x9cSemblance Processing of Borehole Acoustic Array Dataxe2x80x9d, Geophysics, Vol. 49, pp. 274-281 (March 1984). In STC processing, the measured signal is time window xe2x80x9cfilteredxe2x80x9d and stacked, and a semblance function is computed. The semblance function relates the presence or absence of an arrival with a particular slowness and particular arrival time. If the assumed slowness and arrival time do not coincide with that of the measured arrival, the semblance takes on a smaller value. Consequently, arrivals in the received waveforms manifest themselves as local peaks in a plot of semblance versus slowness and arrival time. These peaks are typically found in a peak-finding routine discussed in the aforementioned article by Kimball and Marzetta.
Prior art FIG. 1 illustrates an array of eight waveforms recorded by a sonic tool at a given depth. Such an array is referred to as xe2x80x9ca framexe2x80x9d. Each frame is processed using STC. The first stage in STC is stacking or beamforming of the waveform, to compute semblance, a two-dimensional function of slowness and time which is illustrated in prior art FIG. 2 and is generally referred to as the STC slowness-time plane. The semblance function is given by Equation (1) where xi(t) is the waveform recorded by the i-th receiver of an M-receiver equally spaced array with inter-receiver spacing xcex94z. The array of waveforms {xi(t)} acquired at depth z constitutes a single frame of data.                               ρ          ⁡                      (                          τ              ,              p                        )                          =                                            ∫                              τ                -                                  τ                  2                                                            τ                +                                  τ                  2                                                      ⁢                                                            "LeftBracketingBar"                                                            ∑                                              k                        =                        0                                                                    M                        -                        1                                                              ⁢                                          xe2x80x83                                        ⁢                                                                  x                        i                                            ⁡                                              (                                                  t                          -                                                      k                            ⁢                                                          xe2x80x83                                                        ⁢                            Δ                            ⁢                                                          xe2x80x83                                                        ⁢                            zp                                                                          )                                                                              "RightBracketingBar"                                2                            ⁢                              xe2x80x83                            ⁢                              ⅆ                t                                                          M            ⁢                                          ∫                                  τ                  -                                      τ                    2                                                                    τ                  +                                      τ                    2                                                              ⁢                                                ∑                                      k                    =                    0                                                        M                    -                    1                                                  ⁢                                                                            "LeftBracketingBar"                                              xe2x80x83                                            ⁢                                                                        x                          i                                                ⁡                                                  (                                                      t                            -                                                          k                              ⁢                                                              xe2x80x83                                                            ⁢                              Δ                              ⁢                                                              xe2x80x83                                                            ⁢                              zp                                                                                )                                                                    "RightBracketingBar"                                        2                                    ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                    t                                                                                                          (        1        )            
The semblance xcfx81(xcfx84,p) for a particular depth z is a function of time xcfx84 and slowness p. More particularly, it is the quotient of the beamformed energy output by the array at slowness p (the xe2x80x9ccoherent energyxe2x80x9d) divided by the waveform energy in a time window of length T (the xe2x80x9ctotal energyxe2x80x9d).
The second step in STC processing is to identify peaks in the slowness-time plane with waveform arrivals propagating across the array. Peaks are identified by sweeping the plane with a peak mask. The peak mask is a parallelogram having a slope which corresponds to the transmitter-receiver spacing. A peak is defined as a maximum within the mask region. For each peak, five variables are recorded: the slowness coordinate p, the time coordinate xcfx84, the semblance xcfx81(xcfx84,p), the coherent energy (the numerator of Equation 1), and the total energy (the denominator of Equation 1).
The third step in STC processing is called xe2x80x9clabelingxe2x80x9d or xe2x80x9cclassificationxe2x80x9d and this involves a classification of the peaks found in the second step. The peaks will be classified as corresponding to compressional (P-wave), shear (S-wave) or Stoneley arrivals. Correct classification of the peaks is a difficult process for a number of reasons. Some of the peaks may correspond to spatial aliases rather than the arrival of real waveforms. Some of the peaks may actually be two peaks close together. In general, the problem with state-of-the-art labeling is that small changes in data can cause large differences in the final classification. For example, one of the measures used to classify peaks is the xe2x80x9cslowness projectionxe2x80x9d which is given by Equation (2) and illustrated in prior art FIG. 3. The slowness projection I(p,z) is essentially slowness as a function of depth z which is found at maximum semblance over time xcfx84.                               I          ⁡                      (                          p              ,              z                        )                          =                              max            τ                    ⁢                      ρ            ⁡                          (                              τ                ,                p                ,                z                            )                                                          (        2        )            
The slowness projection (FIG. 3) is currently used as a qualitative check on the slowness logs (prior art FIG. 4) output by the labeling function. Bright tracks in the slowness projection may be identified with the arrival of waveforms. The logs are superimposed on the slowness projections to determine whether the logs are consistent with the slowness projections. In certain regions where the arrival is weak or absent, the intensity of the track diminishes or becomes diffuse. For example, as shown in FIG. 4, the compressional and shear arrivals in a particular borehole seem to vanish at about 3400 feet and deeper. According to actual state-of-the-art algorithms, the classification of shear slowness fails at approximately 3410 feet and the classification of compressional arrival fails at approximately 3425 feet. While abrupt changes in the output of the labelling function may be obvious to the eye, they are difficult to detect automatically.
In addition, the following documents are incorporated herein by reference in their entirety: (1) Harvey, A. C., Time Series Models, MIT Press (1994), Chapter 4 and (2) Scharf, L., Statistical Signal Processing, Wiley (1991), pages 306-311 and 464-469.
It is therefore an object of the invention to provide methods for more accurately labeling sonic waveform information.
It is also an object of the invention to provide methods for tracking sonic measurements into sequences which may be identified as belonging to a single arrival or xe2x80x9ctrackxe2x80x9d.
It is a further object of the invention to provide methods for classifying tracks as representative of a particular type of arrival.
It is another object of the invention to provide methods of waveform analysis which can be performed automatically.
In accord with these objects which will be discussed in detail below, the methods of the present invention include a tracking component and a classification component. More particularly, the methods employ xe2x80x9ctracking algorithmsxe2x80x9d and Bayesian analysis to classify sonic waveforms subjected to STC processing. According to the tracking part of the invention, a probability model is built to distinguish true xe2x80x9carrivalsxe2x80x9d (e.g. compressional, shear, etc.) from xe2x80x9cfalse alarmsxe2x80x9d (e.g. noise) before the arrivals are classified. The probability model maps the xe2x80x9ccontinuityxe2x80x9d of tracks (slowness/time over depth) and is used to determine whether sequences of measurements are sufficiently xe2x80x9ccontinuousxe2x80x9d to be classified as tracks. The probability model is used to evaluate the likelihood of the data in various possible classifications (hypotheses). The classification part of the invention seeks to find the hypothesis with the maximum xe2x80x9cposterior probabilityxe2x80x9d according to Bayes"" Theorem. Prior probabilities and posterior probabilities are constructed for each track using the tracking algorithm and the posterior probabilities are evaluated by the classification part of the method. The posterior probabilities are computed sequentially and recursively, updating the probabilities after each measurement frame at depth k is acquired.
According to a presently preferred embodiment of the tracking algorithm, a data set from each frame includes the measurements paired with an indication of the number of measurements in each frame. The measurements are considered to be xe2x80x9ccontinuousxe2x80x9d whereas the number of measurements is discrete. Separate probability models are constructed for the continuous and discrete data.
The discrete probability model is based upon the xe2x80x9cmodel orderxe2x80x9d O, which is the number of tracks, and the number of measurements implied by the N hypotheses, where N is the sum of the number of tracks observed in the frame and the number of false alarm measurements. Hypotheses are calculated for each frame in a set of k frames where the hypotheses for each frame is a function of the current model order, the previous model order implied by the ancestor hypotheses, and the number of measurements implied by the previous frames. The hypotheses completely determine the possible classifications of measurements in the current frame and their relationship to measurements in the previous frame. The tracks are ranked according to increasing slowness. All of the measurements are classified as to which track they belong and tracks containing no measurements are so classified. A matrix of the classifications is created in which each column corresponds to a track and each row corresponds to a measurement. Probabilities are assigned to each classification based on the model orders and the number of tracks observed as well as the number of false alarms, the number of false alarms being modeled as a Poisson distribution with a known mean. Another matrix is created where the first column specifies which tracks in the previous frame may correspond to track 1 in the current frame, the second column specifies which tracks in the previous frame may correspond to track 2 in the current frame, etc. Probabilities are assigned to the cells of the second matrix based on the classifications, the track assignments, and the probability models described above. A prior probability model for the discrete hypotheses together with the data likelihood are then created based on the foregoing.
The continuous probability models (the likelihood that a sequence of measurements is sufficiently continuous to qualify as a track) are based on the slowness-time coordinates of every peak. A Kalman filter model is used for the slowness-time coordinates. Slownesses are correlated across depth by modeling the sequence of slownesses as the output of an ARMA filter with known coefficients, driven by white Gaussian noise. The slownesses are not assumed to be correlated with times. The actual probability computations are implemented with Kalman filters. Each slowness-time pair from an STC peak is modeled as a two-dimensional vector. A two-dimensional Kalman filter is made from two one-dimensional filters, one for the slowness sequences and one for the time sequences. No correlation between slowness and time is assumed for this part of the method. The Kalman equations enable the evaluation of the likelihood of the sequence of measurements forming a track given the prior continuity model for the track. At each depth, Kalman prediction gives the two-dimensional Gaussian probability distribution of slowness-time data given the past measurements (previous data). The likelihood of the data belonging to the track is computed by evaluating the Gaussian distribution at the slowness-time point specified by the data. A Kalman update of the state and covariance using the data at this depth predicts the two-dimensional probability distribution of the slowness-time at the next depth.
In the classification part of the method, the discrete and the continuous parts of the probability model are subjected to Bayesian analysis. More particularly, the posterior probabilities for each frame are evaluated using Bayes"" Theorem and the most probable hypothesis is selected to classify the tracks. Three different probability models may be used: xe2x80x9cpessimisticxe2x80x9d, xe2x80x9cforwardxe2x80x9d, and xe2x80x9cforward/backwardxe2x80x9d. Pessimistic classification uses only the data in the current frame for classification to compute posterior probability. Forward classification incorporates information provided by all data up to and including the current frame. Forward/backward classification, which is presently preferred, incorporates information from some number of prior frames and some number of subsequent frames. Forward classification can operate in real time with substantially no latency. Forward/backward classification can be run with a latency depending on the number of frames included in each classification. Several assumptions are made to simplify the analyses.
Additional objects and advantages of the invention will become apparent to those skilled in the art upon reference to the detailed description taken in conjunction with the provided figures.